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A recently proposed connection between the threshold for the stability of the iterative 
solution of integral equations for the pair correlation functions of a classical fluid and 
the structural instability of the corresponding real fluid is carefully analyzed. Direct 
calculation of the Lyapunov exponent of the standard iterative solution of HNC and 
PY integral equations for the ID hard rods fluid shows the same behavior observed 
in 3D systems. Since no phase transition is allowed in such ID system, our analysis 
shows that the proposed one phase criterion, at least in this case, fails. We argue 
that the observed proximity between the numerical and the structural instability in 
3D originates from the enhanced structure present in the fluid but, in view of the 
arbitrary dependence on the iteration scheme, it seems uneasy to relate the numerical 
stability analysis to a robust one-phase criterion for predicting a thermodynamic 
phase transition. 

PACS numbers: 05.70.Fh, 61.20.Ne 



I. INTRODUCTION 

When studying the structure and thermodynamics of classical fluids one is often faced 
with the task of solving the nonlinear integral equation which stems out of the combination 
of the Ornstein-Zernike equation and an approximate relation between pair potential and 
correlation functions (the closure) jl|. Integral equations can be generally written in the 
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form 

7 (r) = A 7 (r) , (1) 

where 7 (r) G 5* may be the total correlation function h(r), the direct correlation function 
c(r), or a combination of the two, S is a set of a metric space of functions, and A : S — > S 1 
is a non linear operator mapping 5* into itself. 

Numerical analysis of integral equations suggests the use of the following combination 

7 (r) = h(r) - c(r) , (2) 

since 7 is a much smoother function than h or c, especially in the core region. 

It has been pointed out by Malescio et. al. that, amongst the different numerical 

schemes that one may choose to solve the simple iterative scheme of Picard plays a 
special role. Picard scheme consists in generating successive approximations to the solution 
through the relationship 

7n+l = , (3) 

starting from some initial value 7 o- If the sequence of successive approximations {7 n } con- 
verges toward a value 7 *, then 7 * is a fixed point for the operator A i.e. it is a solution of Eq. 

7* = A'j*. Banach' s fixed point theorem (see chapter 1 in p| especially theorem l.A) 
states that, given an operator A : S — > S, where S is a closed nonempty set in a complete 
metric space, the simple iteration (J3J) may converge toward the only fixed point in S (A is 
k contractive) or it may not converge (A is non expansive). So the simple iterative method 
can be used to signal a fundamental change in the properties of the underlying operator. 

The operator A will in general depend on the thermodynamic state of the fluid. In order 
to determine the properties of the operator at a given state we can proceed as follows. First, 
we find the fixed point 7 * using a numerical scheme (more refined then Picard' s) capable 
of converging in the high density region. Next, we perturb the fixed point with an arbitrary 
initial perturbation S (r) so that 

dA 



A( 7 * + 5 ) ~ A^ + o _ 



5 = Y + M5 , (4) 



where we have introduced the Floquet matrix M. Now 5\ = M5q may be considered as the 
new perturbation. We then generate the succession {5 n } where 

5 n = M5 n _! . (5) 
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If the succession converges to zero then the operator A is k contractive, if it diverges the 
operator is non expansive. Malescio et. al. call {5 n } fictitious dynamics and associate to 
the resulting fate of the initial perturbation the nature of the structural equilibrium of the 
fluid. If the succession converges to zero they say that the fluid is structurally stable and 
structurally unstable otherwise. We will call pi ns t the density where the transition between 
a structurally stable and unstable fluid occurs. 

Following Malescio et. al. it is possible to define a measure for the structural stability of 
the system as follows. We define 

_ ||Mft(r)|| 

1 \Mr)\\ ' [) 



where ||/(r)|| = y YliLi / 2 ( r «) is the norm of a function / defined over a mesh of N points. 
We assume that the norm of the perturbation depends exponentially on the number of 
iterations 

ll<U = IN|2 A " , (7) 

where A is the Lyapunov exponent related to the fictitious dynamics. Then one can write 
the average exponential stretching of initially nearby points as 



A 



lim-log 2 [T\SA . (8) 

\i=0 / 



Malescio et. al. have calculated the dependence of A on the density for various simple 
three dimensional liquids (and various closures): hard spheres Yukawa, inverse power 
and Lennard- Jones potentials |3j. For all these systems they found that A increases with 
the density and the density at which A becomes positive, Pi ns t, falls close to the freezing 
density p/ of the fluid system. This occurrence lead them to propose this kind of analysis 
as a one-phase criterion to predict the freezing transition of a dense fluid and to estimate 
Pf. However, we think that there are some practical and conceptual difficulties with such 
one-phase criterion. 

First of all, it does not depend only on the closure adopted but also on the kind of 
algorithm used to solve the integral equation. Indeed, different algorithms give different 
Pi ns t and Malescio et. al. choose to use as instability threshold for their criterion the one 
obtained using Picard algorithm, thus giving to it a special status. However, it is hard to 
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understand why the particular algorithm adopted in the solution of the integral equation 
should be directly related to a phase boundary. 

Moreover, one would expect that the estimate of pi ns t would improve in connection with 
improved closures. This is not the case, at least in the one component hard sphere fluid. 

Even a more serious doubt about the validity of the proposed criterion comes from its 
behavior in one dimensional systems. In this paper we present the same Lyapunov exponent 
analysis on a system of hard rods in one dimension treated using either the Percus-Yevick 
(PY) or the hypernetted chain (HNC) approximations. What we find is that the Lyapunov 
exponent as a function of density has the same behavior as that for the three dimensional 
system (hard spheres): it becomes positive beyond a certain pi ns t- Since it is known |y] that 
a one dimensional fluid of hard rods does not have a phase transition, our result sheds some 
doubts on the validity of the proposed criterion. 



As numerical scheme to calculate the fixed point we used Zerah' s algorithm [7( for 
the three dimensional systems and a modified iterative method for the hard rods in one 
dimension. In the modified iterative method input and output are mixed at each iteration 



where a is a real parameter < a < 1. Note that while for a non expansive operator A 
the Picard iterative method (jHJ) needs not converge, one can prove convergence results on 
an Hilbert space for the modified iterative method with fixed a (see proposition 10.16 in 
j^). In all the computations we used a uniform grid of N = 1024 points with a spacing 
5r = 0.025. Generally, we observed a marginal increase of pi nst by lowering N. 

A method to find a Lyapunov exponent, equivalent but more accurate than the one of 
Malescio et. al. (jHJ), goes through the diagonalization of the Floquet matrix. Note that 
in general this matrix is non symmetric, thus yielding complex eigenvalues. A Lyapunov 
exponent can then be defined as [s| 



where erj and e% are respectively the real and imaginary part of the i-th eigenvalue. In our 



II. TECHNICAL DETAILS 
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numerical computations we always used recipe (fTUj) to calculate the Lyapunov exponents 
since it is explicitly independent from the choice of an initial perturbation. 

We constructed the Floquet matrix in the following way j2(. In a Picard iteration we 
start from j(r) we calculate c(r) from the closure approximation, we calculate its Fourier 
transform c(k), we calculate j(k) from the OZ equation, and finally we anti transform 7 to 
get 7'(V)- For example for a three dimensional system a PY iteration in discrete form can 
be written as follows 



(1 + 7.) 



1) 



N-1 



i=l 

,s2 



Ti sin(/c J rj)cj 



li = PCj/Q-pZj) 



5k 



N-1 



li 



2ir 2 r 



- ^ k i sin(A; i r i )7 i 



(11) 
(12) 

(13) 

(14) 



where = i5r are the N mesh points in r space, kj = j5k are the N mesh points in k 
space, with 5k = n/(N5r), Ci = c(r,), 7« = 7(r,), Cj = c(fcj), 7, = 7(fcj)? an d 0? = <K r i) is 
the interparticle potential calculated on the grid points. The Floquet matrix will then be 



Mi. 



<H_ = \p Hi dim dc m dc-j 
d lj ~[ ®lm dc m dcj d^/j 
5r5k ( rj 
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(£) (e~^ - 1) (a-, - a + ,; 



(15) 



where 



iV-l 



Di = cos (£7, 



m=l 



2pc r , 



1 - pc r , 



+ 



1 - pc r , 



(16) 



The HNC case can be obtained replacing in ((TSj) [exp(— {3<frj) — 1] with [exp(— (3<j)j+7j) — 1]. 

To derive the expression for the Floquet matrix valid for the one dimensional system and 
consistent with a trapezoidal discretization of the integrals, we need to replace (fT2"|) and (fTljl 
with 



/N-1 1 \ 

Cj = 25r I ^ cos^kjr^Ci + -c j 



li = ~[ ^2 C0S ( k j r i)lj + ^10 



'N-1 



(17) 
(18) 



i=l 
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III. NUMERICAL RESULTS 

We checked our procedure for a three dimensional hard spheres fluid and a Lennard- Jones 
fluid at a reduced temperature T* = 2.74. Our results, obtained using recipe (fTUjl . were in 
good agreement with those of Malescio et. al. which used recipe (jHJ) instead (another 

difference between our analysis and theirs is that we used for 7 the indirect correlation 
function © while they were using the total correlation function h). For the Lennard- Jones 

n 

fluid our results where indistinguishable from those of Malescio et. al. |3J). We found a 
reduced instability density p* nst around 1.09 in the PY approximation and around 1.06 in 
the HNC approximation. For the three dimensional hard sphere fluid we found slightly larger 
(4%) values for pi ns t- We found a r] inst = p inst iTd 3 /6 of about 0.445 in the PY approximation 
and around 0.461 in the HNC approximation. We also checked the value corresponding to 
the Martynov-Sarkisov (MS) jl0| closure and we found f] inst = 0.543. 

We feel that the differences are within what we can expect on the basis of small numerical 
differences in different implementations. We think that it is more worth of notice that 
closures providing better structural and thermodynamic properties, like PY or MS do not 
provide a better value of r] inst . 

We looked at the structure of the Floquet matrix too but from direct inspection we can 
conclude that it is not diagonally dominated. 

Then, we have calculated the Lyapunov exponent ()10j) as a function of the density for 
a fluid of hard rods in one dimension using both PY and HNC closures. The results of 
the calculation are shown in Fig. Q and Fig. El The curves show the same qualitative 
behavior as the ones for the three dimensional fluid. From Fig. ^we can see how the slope 
of the curves starts high at low densities and decreases rapidly with p. At high densities the 
Lyapunov exponent becomes zero at Pi ns t- As expected, to find the fixed point for p > p inst 
it is necessary to choose a < 1 in the modified iterative scheme 0. Before reaching the 
instability threshold the curves show a rapid change in their slope at p c < pi ns t- Figure 
12] shows a magnification of the region around p c from which we are lead to conclude that, 
within the numerical accuracy of the calculations, the slope of the curves d\' /dp undergoes 
a jump at p c . 
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IV. CONCLUSIONS 

The fictitious dynamics associated to the iterative solution of an integral equation can 
signal the transition of the map of the integral equation from k contractive to non expansive. 
If the Lyapunov exponent is negative the map is k contractive, if it is positive the map is 
non expansive. 

Since it is possible to modify in an arbitrary way the fictitious dynamics keeping the 
same fixed point, it is difficult to understand a deep direct connection between the stability 
properties of the map and a one-phase criterion for a thermodynamic transition. 

Admittedly the correlations shown by Malescio et al. are striking. We calculated the 
Lyapunov exponent as a function of the density for various fluids (hard spheres in one 
and three dimensions and three dimensional Lennard- Jones fluid) both in the HNC and 
PY approximations. For the three dimensional fluids the instability density falls close to 
the freezing density pf. For example, the Lennard- Jones fluid studied with HNC should 
undergo a freezing transition at p* ~ 1.06 or at p* ~ 1.09, if studied with PY , rather 
close to the freezing density pj- ~ 1.113. For hard spheres p* nst is about 10% smaller than 
Pj ~ 0.948. The Hansen- Verlet "rule" states that a simple fluid freezes when the maximum 
of the structure factor is about 2.85 [llj. According to this rule the three dimensional 
hard spheres fluid studied with HNC should undergo a freezing transition at p ~ 1.01 while 
when studied with PY the transition should be at p ~ 0.936. The corresponding estimates 
obtained through p* nst , 0.879 (HNC) and 0.850 (PY) are poorer and, more important, are 
not consistent with the well known better performance of PY in the case of hard spheres. 
In one dimension, a fluid of hard spheres (hard rods), cannot undergo a phase transition 
From Fig. ^ we see that the system still becomes structurally unstable. This can 
be explained by observing that the structural stability as defined by Malescio et. al. is 
a property of the map A and in particular of the algorithm used to get solution of the 
integral equation under study. As such, it is not directly related to the thermodynamic 
properties even at the approximate level of the theory (there is no direct relation between 
the contractiveness properties of A and the thermodynamics). It looks more reasonable that 
the increase of the correlations would be the common origin of the numerical instability of 
the Picard iteration and, whenever it is possible, of thermodynamic phase transitions. 
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LIST OF FIGURES 

1 For a fluid of hard rods in one dimension, we show the Lyapunov exponent as a function 
of the reduced density (p* = pa where a is the rods width) as calculated using the PY 
and the HNC closures. 

2 We show a magnification of Fig. ^ in a neighborhood of the instability threshold. 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



P 

FIG. 1: R. Fantoni and G. Pastore 




-0.10 1 ■ ■ ■ ■ 1 

0.60 0.62 0.64 0.66 0.68 0.70 

* 

P 



FIG. 2: R. Fantoni and G. Pastore 



